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Abstract. 

The effects of the ion Larmor radius on magnetic reconnection are investigated by means 
of numerical simulations, with a Hamiltonian gyrofluid model. In the linear regime, it is found 
that ion diamagnetic effects decrease the growth rate of the dominant mode. Increasing ion 
temperature tends to make the magnetic islands propagate in the ion diamagnetic drift direction. 
In the nonlinear regime, diamagnetic effects reduce the final width of the island. Unlike the 
electron density, the guiding center density does not tend to distribute along separatrices and 
at high ion temperature, the electrostatic potential exhibits the superposition of a small scale 
structure, related to the electron density, and a large scale structure, related to the ion guiding- 
center density. 



1. Introduction 

Magnetic reconnection is responsible for many phenomena occurring in laboratory and 
astrophysical plasmas. [TJ [2] In this paper, we restrict consideration to reconnection in low- 
beta configurations (i.e. m e /Mj < j3 < 1, with m e and Mi indicating electron and ion mass, 
respectively ) characterized by the presence of a magnetic "guide" field perpendicular to the 
reconnection plane. The guide field ensures that the ions undergo rapid Larmor gyration as 
they cross the reconnection region. Guide field reconnection is thought to be responsible, in 
particular, for the sawtooth crash in tokamaks. Kinetic simulations of guide-field reconnection 
are generally more costly and difficult than in the alternative case of "component reconnection," 
due to the separation of scales imposed by the guide field. [31 HI [5] In the linear [6] as well as 
in the nonlinear regime, [HJ [9] magnetic reconnection is characterized by processes occurring 
in very narrow resonant layers. For high-temperature plasmas, such as those existing in the 
present generation of tokamaks, the resonant layer is narrower than the ion Larmor radius, so 
that the ions may traverse the entire layer during a cyclotron period. [10] Their response is then 
a nonlocal functional of the fields they sample along their orbit. |11[ [12] 

Conventional fluid models treat finite Larmor radius (FLR) effects by expanding the ion 
response to first order in (k±pi) 2 . [131 E2 ESI US] In the linear regime, the predictions of fluid 
models that include FLR effects [IT] are qualitatively similar to the results of kinetic calculations 
|11[ [12] that account for the nonlocal response of the ions. In particular, both fluid and kinetic 



models predict the stabilization of ideal instabilities when the growth rate 70 computed in 
the absence of diamagnetic drifts is less than half of u;*j, the ion diamagnetic drift angular 
frequency. [11[ \T7\ [TB] It is unclear a priori, however, whether fluid FLR calculations such 
as the investigations of the nonlinear saturation of reconnecting modes in Refs. [13|. \15\ I19j . 
offer an adequate description of ion dynamics in the nonlinear regime. In particular, it is 
unclear whether they are adequate for the task of extrapolating the experimental observations 
to different values of pt/L, where L describes the plasma size. Some indications are provided 
by nonlinear simulations using gyrokinetic codes, [20] but numerical limitations make scaling 
studies with such codes impractical. 

The need for a simpler nonlinear model for the nonlocal ion response has motivated the 
derivation and investigation of gyrofluid models that account for these physical effects in a fluid 
context. Grasso et al. [2T] introduced the first electromagnetic gyrofluid model and applied it to 
the study of magnetic reconnection. Their gyrofluid model, although unable to describe kinetic 
phenomena such as wave particle resonance (Landau damping) , [12] nevertheless illuminated the 
qualitative properties of the reconnection dynamics. More elaborate models that do account for 
Landau damping were subsequently developed by Snyder and Hammett, but the applications of 
their models have been limited to the study of electromagnetic effects on turbulence driven by 
pressure gradients. [22] 

In the present contribution, we investigate a recently derived 3-field gyrofluid model |23] by 
means of numerical simulations. This gyrofluid model is part of a family of Hamiltonian fluid 
models [Ml Ell ESI [26] that have been investigated in past years, but is distinguished from 
these other Hamiltonian models by its nonlocal treatment of the ion density evolution. In the 
linear regime, the nonlocal description of the ion density evolution improves on FLR models 
by describing the evanescence of drift waves excited at frequencies in the ion drift direction 
but below the diamagnetic frequency. [23] The propagation and evanescence properties of drift 
waves can be important for determining the role of the polarization current in island evolution. 



In recent work, [28] we have used the 3-field model of Ref. [23J to revisit earlier investigations 
|21| of gyrofluid magnetic reconnection in a homogeneous plasma, where ion diamagnetic effects 
vanish. The present article extends these studies to the case of inhomogeneous plasma, where ion 
diamagnetic drifts affect the reconnection. Such ion drift effects can be important for sawtooth 
dynamics in tokamaks \13\ 115] . where it accounts for the saturation of the resistive kink 
instability which can prevent complete reconnection during sawtooth crashes. \29\ [30] It is also 
important for reconnection at the dayside magnetopause [31[ 132] where the interface between 
the incoming solar wind and the Earth's magnetosphere is characterized by a strong density 
gradient. Unfortunately, our assumption of a strong guide field makes our model inapplicable 
to the magnetopause. 

The paper is organized as follows. In Sec. [2] we review the model equations. Sec. is devoted 
to the simulation results in the linear phase and to the comparison with analytical predictions. 
In Sec. H] we consider the evolution of the system in the nonlinear phase. We conclude in Sec. EJ 

2. The model 

We consider the two-dimensional version of the Hamiltonian gyrofluid model described in 
Ref. [53] . The model equations are 
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n e = rl /2 m + (r - i^/pf. (4) 

These represent the continuity equation for the ion guiding centers, the electron continuity 
equation, the Ohm's law and the Poisson's equation, respectively. In Eqs.([T|)-(j3]) rn is the ion 
guiding center density, n e the electron density, 4> the electrostatic potential, i[) the poloidal 

1/2 

magnetic flux function, d e the electron skin depth, p s the sonic Larmor radius and = r o (/> 
is the gyro-averaged electrostatic potential. For our analysis the Pade approximant version 
Tl /2 = (l-pfV 2 /2)- 1 of the gyro-average operator will be adopted, where pi is the ion Larmor 
radius. Given a Cartesian coordinate system (x, y, z), we assume all the fields be translationally 
invariant along z and we define the canonical bracket between two generic fields / and g by 
[/, g] = z ■ V/ x V<?. We also recall that we normalize the time with respect to the Alfven time 
and distances with respect to a magnetic equilibrium scale length L. Dependent variables are 
normalized in the following way 

L hi Lh e A z p 2 s L ecj) 

di no di n BL L l d { T e 

where carets denote dimensional variables, di is the ion skin depth, uq a background density 
amplitude, A z the magnetic potential, e the unit charge, B a characteristic magnetic field 
amplitude and T e the electron temperature, which is assumed to be constant. 



3. Numerical simulations in the linear regime 

We devote this section to the analysis of the numerical solutions of the model {!])-© i n the 
phase in which nonlinear terms are not dominant. 

We solve the model equations over the domain {(x,y) : —ir < x < ir, —air < y < air}, where we 
prescribe the value of the number a at each simulation. The grid is composed by 1024 x 128 points 
and double periodic boundary conditions are imposed on the field perturbations, which are the 
quantities that the code advances in time. Unlike Ref. [28J we consider now inhomogeneous 
density equilibria. The simulations are then started by perturbing the equilibrium 

11 

%W = V- n £eq (x) = n' x, ip eq {x)= ^ a n exp(inx), (5) 

n=-ll 

where n' is a constant and the a n are the Fourier coefficients of the function f(x) = 1/ cosh 2 x. 
With respect to the simulations of Ref. [25] the presence of equilibrium density gradients 
makes the evolution of rij no longer negligible and introduces the diamagnetic effects. A global 
dispersion relation for the system in the presence of diamagnetic effects was derived in Ref. |23j . 
This showed also how the Pade approximant version of the gyrofluid model was able to reproduce 
adequately dispersive properties of the system, such as a spectral gap in the frequencies, which 
are not caught by Finite Larmor Radius (FLR) models. 

The equilibrium ([5]) is perturbed in raj with a four-cell pattern disturbance of the form 
hi (x cos(x + y/a) — cos{x — y/a). The field <fi is also perturbed according to (jU), in such a 
way that the initial perturbation on n e is zero. 

If we indicate with xi. x iVi^) an y °f the three dynamical variables of our system, i.e. raj, 
n e and ip, then we can decompose the fields as \ = Xeq(x) + x(x,y,t). We write then the 
perturbation \ as a Fourier series truncated at N modes, so that \ = J2-N Xk ex P(^ ' x + 
where the hat indicates the Fourier transform, k is the wave vector and A = 7 — ioj is a complex 
number having the growth rate 7 as real part and the rotation frequency — u) as imaginary 
part. Simulations with a relatively small amplitude of the initial perturbation, corresponding 
to 10 -6 , have been carried out, in order to obtain a long linear phase (~ 25 Alfven times for 



Table 1. Table displaying values of linear growth rate 7 and rotation frequency u, for the 
dominant mode. A comparison is made between values obtained from numerical simulations 
(%™,7™m) and values (u} theor ,^ theor ) predicted by the asymptotic theory of Ref. p]. The 
agreement between numerical and analytical results is better for A' = 59.9, which is a value 
closer to the asymptotic regime of validity of the theory. For all cases d e = 0.2. 
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most of the simulations) that can be effectively investigated. In particular, we are interested 
in extracting from the simulations the values of 7 and cj, which can permit us to observe the 
effects of equilibrium density gradients and finite gyro-radius, on the linear mode stabilization 
and island rotation frequency. We find it also of interest to compare the obtained numerical 
results with those predicted by the theory developed in Ref. [11]. The linear analysis by Porcelli 
refers namely to the collisionless case in the presence of finite ion temperature. Here we limit 
ourselves to recall the results of the theory that are relevant for our analysis. By adopting 
a matching asymptotic technique, in Ref. [11], the following dispersion relation, for the main 
mode, is obtained 

A 2 - w* e w*j + i\(uj*e + « 7q. (6) 

In §§§ the electron diamagnetic frequency u>* e is defined by u>* e = — k y cT e n' / (eBnoVA), where va 
is the Alfven speed based on B. The ion diamagnetic frequency is u*i = —ruj* e , with r = pf/p 2 , 
and 70, for our case, reads 70 = 2k y (2d e (p 2 s + p 2 )/-7r) 1//3 . Splitting the dispersion relation into its 
real and imaginary part yields 

«»^(1-t), (7) 
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which give us a prediction for the growth rate and the rotation frequency. Such results, however, 
are valid only in the limit pi 3> d e and A' ^> p7 d e , where A' is the standard stability 
parameter for tearing modes. 

Examples of theoretical and numerical values are provided in Table [TJ Two main cases are 
considered, corresponding to two different values of A' (which correspond to the values a = 1.5 
and a = 2 of the parameter scaling the length of the domain along y). For each of these two cases 
the same sets of parameters pi, p s and v* = u)* e /k y are considered. From these representative 
examples we can see that the numerical results follow the trends predicted by the theory, for 
the dominant mode (the latter is k y = 0.66 for A' = 34.2 and k y = 0.5 for A' = 59.9). In 
particular, for fixed pi and p s , the cases with v* = —0.4 exhibit smaller growth rate than 
those with v* = —0.16. This indicates the stabilizing effect operated by diamagnetic flows and 
corresponding to the term proportional to w 2 e in ([8]). We report also that another simulation, 
again with pi = p s = 0.4, but with v* = —0.64, showed no growth of the magnetic island, 



indicating the complete linear stabilization due to diamagnetic effects. For these two sets of 
parameters, the numerical solutions show almost no rotation frequency. On the other hand, for 
the case pi = 0.35, p s = 0.2 and v* = —0.04, the numerical solutions show a positive value of u, 
which corresponds to a propagation along the direction of the ion diamagnetic drift (which in 
our case corresponds to the positive y direction). The propagation along the ion diamagnetic 
drift direction is indeed also what is predicted by the theory, according to ([7]), when ions are 
hotter than electrons. With regard to the quantitative agreement with the theory, Table[T]shows 
that this improves when increasing A'. In the case with larger A' the distances between the 
numerical and predicted values get shorter and we can quantify that the relative error drops 
from an average of 28% to 18.7%. The reduction of the error is of course a consequence of the 
fact that, by increasing A', we get closer to the asymptotic regime in which the theory is valid. 

4. Numerical simulations in the nonlinear regime 

In this section we focus on the evolution of the system in the nonlinear phase. 
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Figure 1. Contour plots of tft for v* = —0.16, at t = 39 (left) and = —0.4 at t = 50 (right). 
The values of the other parameters are: p s = pi = 0.4, d e = 0.2. The plots indicate that stronger 
diamagnetic effects lead to a smaller saturated island. 

Fig. [H shows a comparison between contour plots of tp for = —0.16 and = —0.4, at 
times at which the two simulations reached the same phase of the dynamical evolution. From 
the figure we can see that the presence of stronger diamagnetic effects, apart from reducing the 
linear growth rate, also leads to a nonlinear island of smaller width. We recall that a nonlinear 
stabilizing effect of n' was observed in a cold- ion, three-field model in cylindrical geometry |15j . 
In that context diamagnetic flows were indeed suggested to be responsible for the saturation 
of islands at a small width. We also note that, in spite of the fact that, in the linear phase, 
the islands exhibited very little rotation for these parameters, in the nonlinear regime, upward 
shifts of the O-points are present, in particular in the case v * = —0.4. Island propagation in the 
ion diamagnetic drift direction, caused either by viscosity [6'6\ [TB] or by nonlinear zonal flows, 
[Ml ES] has been reported in FLR simulations of visco-resistive reconnection. 
An effective way to look at the nonlinear dynamics of our system is to express the model in 
terms of its "normal fields" , suggested by the Hamiltonian structure. We recall that the model 
equations (HI)-© can be rewritten as [23] 

5? + [*,*]=<), (9) 



-^± + tt + , G+ } = 0, (10) 
r>C 

-^ + [0_,G_]=O, (11) 

where G± = ip — d^V 2 ^ ± d e p s n e and (j>± = 4> ± {Ps/d e )tp. The fields and G±, which we 
denote as normal fields, represent natural variables for the system because, in terms of them, 
the model equations take the form of simple advection equations. As pointed out in a number 
of previous papers on reconnection based on Hamiltonian fluid models [361 EH GS GSl [261 HO] , 
the stream functions <j>± correspond to velocity fields that rotate in opposite directions. The 
fields G±, adverted by such velocity fields, get then stirred in opposite directions. The stirring 
process causes G± to form finer and finer structures, leading to an energy cascade toward small 
scales. Because n e is proportional to the difference between G+ and G_ , such cascade inevitably 
reflects on the electron density structures. 

In order to investigate the nonlinear structures, we double the number of grid points in the 
y direction, which enables us to better resolve small scale structures that form nonlinearly. 
We also bring up to 10~ 3 the amplitude of the initial perturbation. From Fig. [2] indeed we 
see that electron density, as is typical of Hamiltonian reconnection, develops lobes around 
the separatrices. These emerge as the difference between the spiral arms present in G+ and 
G- (the latter not shown here but such that G-(x,y,t) = G+(—x,y,t)). Unlike the above 
mentioned studies of Hamiltonian reconnection, however, the structures of Fig. [2] are of course 
not symmetric with respect to the y = axis, because of diamagnetic effects. The pattern 
of (f) + (and, analogously of 4>-(x,y,t) = — (f>+(— x, y, t)) are such to create thinner and shorter 
layers in the lower part of the magnetic island whereas wider and longer lobes are present in the 
upper part of the magnetic island. Inside the island, stirring of G± takes place and the initial 
linear density profile gets flattened. The ion guiding center density gradient also gets partly 
flattened in the center of the island, as a consequence of the advection operated by <F Such 
stream function, however, has a circulation pattern quite different from that obtained by the 
combination of cj> + and 0_. Therefore does not tend to accumulate around the separatrices, 
unlike n e . The four convective cells of the gyroaveraged electrostatic potential, on the other hand, 
tend to lie along the magnetic separatrices. Such behavior is also evident in the electrostatic 
potential <f). By virtue of the relation <fi = (1 — /9?V 2 /2)$, we can also infer that <fi exhibits, in 
general, finer structures, compared to <5. We can indeed observe regions with steep gradients 
of electrostatic potential, around the separatrices, similarly to what happens to the electron 
density. Nevertheless, in this respect, it is important to recall that, for n' = 0, the growth of 
rii was negligible, and one had 0k ~ —pt^ek-, as pi — > +oo. This led to conclude [2H] that the 
filamentation in n e implied an analogous phenomenon in 0, with consequent formation of strong 
electric fields. In the presence of diamagnetic effects, on the other hand, Eq. (£[]), in the large pi 
limit, yields 

2 

<Ak ~ -P?n ek + ^2^ik- (12) 

Therefore, <j) is no longer proportional to n e , but is given by the superposition of two 
contributions: one, related to h e , which is responsible for the formation of steep gradients, 
and a second one, due to — 2V~ 2 fi; which, on the contrary, is relevant at large scales. 



5. Conclusions 

Numerical simulations of a Hamiltonian electromagnetic gyrofluid model have been carried out 
and investigated. The present analysis represents a step forward, with respect to a recent related 
work [28J, because it accounts for inhomogeneous density equilibria, which lead to diamagnetic 
effects. In the linear phase, diamagnetic effects are shown to stabilize the growth rate of the 




Figure 2. Contour plots of n e , cp, Hi, G + and (f> + at t = 19, for v* = —0.16, Pi = p s = 0.4 
and d e = 0.2. The magnetic island at t = 19 is superimposed onto each plot. 

dominant mode and induce a propagation of the magnetic island. Simulations also indicate 
that the direction of propagation tends to move from that of the electron to that of the ion 
diamagnetic drift, as the ion temperature is increased. Such features are in agreement with what 
is predicted by the linear theory of Ref.|ll] and the quantitative agreement improves as the values 



of parameters approach the asymptotic regime of validity of the theory. In the nonlinear phase, 
stronger equilibrium gradients tend to yield smaller "saturated" islands, as already observed in 
Ref.|15|. We note also a different distribution of the ion guiding center density with respect to 
the electron density, with the former retaining a four-cell pattern, whereas the latter concentrates 
around separatrices, although in an asymmetric way, also with respect to the y = axis, as a 
consequence of the diamagnetic terms. Compared to the homogeneous equilibrium case [28], the 
potential possesses now a smoother contribution, due to the finite n« fluctuations. We finally 
indicate the derivation of a refined version of the linear theory, and the extension of the model 
to include parallel dynamics, as subjects of work in progress. 
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